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Dual phospho/dephosphorylation cycles, as well as covalent enzymatic-catalyzed modifications 
of substrates are widely diffused within cellular systems and are crucial for the control of complex 
responses such as learning, memory and cellular fate determination. Despite the large body of deter- 
ministic studies and the increasing work aimed at elucidating the effect of noise in such systems, some 
aspects remain unclear. Here we study the stationary distribution provided by the two-dimensional 
Chemical Master Equation for a well-known model of a two step phospho/dephosphorylation cycle 
using the quasi-steady state approximation of enzymatic kinetics. Our aim is to analyze the role of 
fluctuations and the molecules distribution properties in the transition to a bistable regime. When 
detailed balance conditions are satisfied it is possible to compute equilibrium distributions in a closed 
and explicit form. When detailed balance is not satisfied, the stationary non-equilibrium state is 
strongly influenced by the chemical fluxes. In the last case, we show how the external field derived 
from the generation and recombination transition rates, can be decomposed by the Helmholtz theo- 
rem, into a conservative and a rotational (irreversible) part. Moreover, this decomposition, allows to 
compute the stationary distribution via a perturbative approach. For a finite number of molecules 
there exists diffusion dynamics in a macroscopic region of the state space where a relevant transition 
rate between the two critical points is observed. Further, the stationary distribution function can 
be approximated by the solution of a Fokker-Planck equation. We illustrate the theoretical results 
using several numerical simulations. 



I. INTRODUCTION 

One of the most important aspects of biological sys- 
tems is their capacity to learn and memorize patterns 
and to adapt themselves to exogenous and endogenous 
stimuli by tuning signal transduction pathways activ- 
ity. The mechanistic description of this behavior is typ- 
ically depicted as a "switch" that can drive the cell fate 
to different stable states characterized by some observ- 
ables such as levels of proteins, messengers, organelles 
or phenotypes |21j . The biochemical machinery of signal 
transduction pathways is largely based on enzymatic re- 
actions, whose average kinetic can be described within 
the framework of chemical kinetics and enzyme reac- 
tions as pioneered by Michaelis and Menten [tJJ [T7] . The 
steady state velocity equation accounts for the majority 
of known enzymatic reactions, and can be adjusted to 
the description of regulatory properties such as cooper- 
ativity, allostericity and activation/inhibition[10 . Theo- 
retical interest in enzymatic reactions has never stopped 
since Michaelis-Menten's work and has lead to new dis- 
coveries such as zero-order ultrasensitivity [TJ 0] . Among 
various enzymatic processes, a wide and important class 
comprises the reversible addition and removal of phos- 
phoric groups via phosphorylation and dephosphoryla- 
tion reactions catalyzed by kinases and phosphatases. 
The phospho/dephosphorylation cycle (PdPC) is a re- 
versible post-translational substrate modification that is 
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central to cellular signalling regulation and can play a 
key role in the switch phenomenon for several biolog- 
ical processes ([3 [15]). Dual PdPC's are classified as 
homogeneous and heterogeneous based on the number 
of different kinases and phosphatases ([B]); the homo- 
geneous has one kinase and one phosphatase, while the 
heterogeneous has two kinases and two phosphatases. 
Variants of homogeneous and heterogeneous dual PdPC's 
may only have a non-specific phosphatase and two spe- 
cific kinases [3] or, symmetrically, a non-specific kinase 
and two specific phosphatases. The PdPC with the non- 
specific phosphatase controls the phosphorylation state of 
AMPA receptors that mediates induction of Long Term 
Potentiation (LTP) and Long Term Depression (LTD) 
in vitro and in vivo [21 [3j [20]. Recently, several au- 
thors have reported bistability in homogeneus pPdPC 
[SI HI]) as well as those with a non-specific phosphatase 
and two different kinases [2|. The bistable behaviour 
of the homogeneous system is explained on the basis of 
a competition between the substrates for the enzymes. 
The majority of studies on biophysical analysis of phos- 
pho/dephosphorylation cycle have been performed in a 
averaged, deterministic framework based on Michaclis- 
Menten (MM) approach, using the steady state approxi- 
mation. However, recently some authors [5] have pointed 
to the role of fluctuations in the dynamics of biochemical 
reactions. Indeed, in a single cell, the concentration of 
molecules (substrates and enzymes) can be low, and thus 
it is necessary to study the PdPC cycle within a stochas- 
tic framework. A "natural" way to cope with this prob- 
lem is the so-called Chemical Master Equation (CME) 
approach |18) . that realizes in an exact way the proba- 
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bilistic dynamics of a finite number of molecules, and re- 
covers the chemical kinetics of the Law of Mass Action, 
which yields the continuous Michaelis-Menten equation 
in the thermodynamic limit (N — > oo,) using the mean 
field approximation. In this paper we study a stochastic 
formulation of enzymatic cycles that has been extensively 
considered by several authors ( [SJ UJ] ) . The determinis- 
tic descriptions of these models characterize the stability 
of fixed points and give a geometrical interpretation of 
the observed steady states, as the intersection of conic 
curves([12 ). The stochastic description can in fact pro- 
vide further information on the relative stability of the 
different steady states in terms of a stationary distribu- 
tion. We propose a perturbative approach for computing 
the stationary solution out of the thermodynamics equi- 
librium. We also point out the role of currents in the 
transition from a mono-modal distribution to a bimodal 
distribution; this corresponds to bifurcation in the deter- 
ministic approach. The possibility that chemical fluxes 
control the distribution shape suggests a generic mecha- 
nism used by biochemical systems out of thermodynamic 
equilibrium to obtain a plastic behavior. Moreover, we 
show that at the bimodal transition there exists a dif- 
fusion region in the configuration space where a Fokker- 
Planck equation can be introduced to approximate the 
stationary solution. Analogous models have been previ- 
ously studied [H ED 03] for single-step PdPC. 



II. DUAL 

PHOSPORYLATION/DEPHOSPHORYLATION 
ENZYMATIC CYCLES 



The process shown in Fig. ([I]) is a two-step chain of 
addition/removal reactions of chemical groups and may, 
in general, model important biological processes such 
as phenotype switching (ultrasensitivity) and chroma- 
tine modification by histone acetylation/deacetylation as 
well as phospho/dephosphorylation reactions. Without 
loss of generality, we perform a detailed study of the ho- 
mogeneous phospho/dephosphorylation two-step cycles 
(PdPC cycles) where two enzymes drive phosphoryla- 
tion and dephosphorylation respectively. Thus, there is a 
competition between the two cycles for the advancement 
of the respective reactions. 



The deterministic Michaelis-Menten (MM) equations 
of the scheme ([lj with the quasi-steady-state hypothesis 
reads: 




FIG. 1. Scheme of the double enzymatic cycle of ad- 
dition/removal reactions of chemical groups via Michelis- 
Menten kinetic equations as shown in eq in the case of 
phosphoric groups. 
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where A and Ap are the concentrations of the non- 
phosphorylated and double phosphorylated substrates, 
k^ii denote the MM constants and are the maxi- 
mal reaction velocities (i — 1, • • • 4). Let n\ and n 2 de- 
note the molecules number of the substrates A and Ap 
respectively, the corresponding CME for the probability 
distribution p(n 1 ,n 2 ,i) is written 



dp 

dt 



= 9i(ni - l,7i 2 )p(ni - l,n 2 ,t) - g 1 (n 1 ,n 2 )p(n 1 ,n 2 ,t) 

+ n(m + l,n 2 )p(n 1 + l,n 2 ,t) - n(ni, n 2 )p(ni, n 2) t) 
+ g 2 (ni,n 2 - l)p(n 1 ,n 2 - l,t) - ,g 2 (ni, n 2 )p(ni, n 2 , t) 
+ r 2 (ni,n 2 + l)p(m,n 2 + l,t) - r 2 (n u n 2 )p{m,n 2 ,t) 

(2) 

where gj{ni,n 2 ) and rj(rii,n 2 ) are the generation and 
recombination terms respectively, defined as : 
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Nt is the total number of molecules, and we have intro- 
duced scaled constants Km = NrkM- The biochemical 
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meaning is that the enzyme quantities should scale as the 
total number of molecule Nt, to have a finite thermody- 
namic limit Nt — > oo, in the transition rates ([3]). As it 
is known from the theory of one-step Markov processes, 
the CME ([2| has a unique stationary solution p s (rii, n 2 ) 
that describes the statistical properties of the system on 
a long time scale. The CME recovers the Mass Action- 
based MM equation in the thermodynamic limit when 
the average field theory approach applies. Indeed it can 
be shown that the critical points of the stationary distri- 
bution for the CME can be approximately computed by 
the conditions (cfr. eq. (15)) 



ffi (711,712) = ri(m + l,n 2 ) 52(711, n 2 ) = r 2 (711,712 + 1) 

(4) 

whose solutions tend to the equilibrium points of the MM 
equation when the fluctuation effects are reduced in the 
thermodynamic limit as 0(l/yN). As a consequence, 
one would expect that the probability distribution be- 
comes singular, being concentrated at the fixed stable 
points of the equations ([I]), and that the transition rate 
among the stability regions of attractive points is negli- 
gible. However, when a phase transition occurs due to 
the bifurcation of the stable solution, fluctuations are 
relevant even for large Nt and the CME approach is 
necessary. In the next section we discuss the stationary 
distribution properties for the CME 

III. THE STATIONARY DISTRIBUTION 

The stationary solution p s (i%i, 712) of the CME (|2) can 
be characterized by a discrete version of the zero diver- 
gence condition for the current vector J components (see 
Appendix) 

J[ = 9i(ni - l,n 2 )p s (ni - l,n 2 ) - ri{ni,n 2 )p s {m,n 2 ) 
J| = g 2 (n\,n 2 - l)p s (n 1 ,n 2 - 1) - r 2 (ni, n 2 )p s (n 1 , n 2 ) 

(5) 

and the CME r.h.s. reads 

D 1 J s 1 {n ll n 2 ) + D 2 J s 2 {n 1 ,n 2 )=Q (6) 

where we have introduced the difference operators 

Dif(n 1 ,n 2 ) = f(ni + 1, n 2 ) - f(m, n 2 ) 
D 2 f(n 1 ,n 2 ) = f{ni,n 2 + 1) - f(m, n 2 ) (7) 

Due to the commutative properties of the difference op- 
erators, the zero-divergence condition for the current 
is equivalent to the existence of a current potential 
A(ni, n 2 ) such that 

Jl(n\,n 2 ) = D 2 A(m,n 2 ) i%i > 1 
J 2 {n\,n 2 ) = —DiA(m,n 2 ) n 2 > 1 

(8) 

We remark that the r.h.s. of eq. ^ is a discrete version 
of the curl operator. The potential difference A(n' 1: n' 2 ) — 



A(rii,n 2 ) defines the chemical transport across any line 
connecting the states (ni,n 2 ) and (n[,n' 2 ). At the sta- 
tionary state the net transport across any closed path is 
zero and we have no current source in the network. As 
discussed in [THl EH] we distinguish two cases: when the 
potential A(ni,n 2 ) is constant (the so called "detailed 
balance condition" ) and the converse case corresponding 
to a non-equilibrium stationary state. In this case the 
stationary solution p s is characterized by the condition 
•h = J 2 = over all the states, whereas in the other 
case we have macroscopic chemical fluxes in the system. 
When detailed balance holds, simple algebraic manipu- 
lations (see Appendix) result in the following conditions 
for the stationary solutions 

Di In. p a (711,712) = In. ai (711,712) (9) 
D 2 In Ps (711,712) = lna 2 (ni,n 2 ) m + n 2 < N T 

(10) 

where the discrete drift vector field components a* are 
defined by: 

91(711,712) 92(711,712) 

Oi(ni,7l2) = — ~, — r a2(7H,7l2)- 



ri (ni + l,7l 2 ) 



7-2(ni,n 2 + 1) ' 
(11) 



Equations (9) and ( 10 ) imply the existence of a potential 



V(ni, n 2 ) such that 

lnai (711,712) = -DiV(ni,n 2 ) 
\na 2 {m,n 2 ) = -D 2 V{nx,n 2 ) (12) 

Using definition ^ it is possible to explicitly compute 
a set of parameter values for the PdPC cycle, that sat- 



isfy the detailed balance condition ( 12 ) according to the 
relations 
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where the reaction velocities Vm are arbitrary. The sta- 
tionary distribution is given by the Maxwell-Boltzmann 
distribution 



p s (ni,n 2 ) = exp(-V(rai,7i2)) 



(14) 



where the potential V(ni,n 2 ) is computed by integrat- 



ing equation (12) and choosing the initial value V(0,0) 
to normalize the distribution (14). Using a thermody- 



namical analogy, we can interpret the potential differ- 
ence V(ni,n 2 ) — V(0,Q) as the chemical energy needed 
to reach the state (ni, n 2 ) from the initial state (0, 0). As 
a consequence, the vector field (11) represents the work 
for one-step transition along the n\ or n 2 direction. Def- 



inition (14) also implies that the critical points of the 



stationary distribution are characterized by the condi- 
tions 



9\{n\,n 2 ) 



g 2 {ni,n 2 ) 



ri(m + 1,712) 7-2(711, 7i2 + 1) 



= 1 



(15) 



and coincides with the critical points of the MM equa- 
tions. In figure [2] we plot the stationary distributions 
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FIG. 2. Stationary distributions for the A and A P states in the double phosphorylation cycle when detailed balance ( 13 1 
holds with Km% = Km* = 1 and Km 2 = K M3 = 2. In the top figure we set the reaction velocities v Ml = »m 2 = 1 and 
vm 2 — vm 3 = 1-05 (symmetric case), whereas in the bottom figure we increase the vm 2 an d vm 3 value to 1.15. The number of 
molecules is Nt = 40. The transition from a unimodal distribution to a bimodal distribution is clearly visible. 
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40 with K Ml = K Mi = 1 
we consider two symmetric 
= 1.05 or 
In 



in the case Nt 
ana Km 2 = K Ms = 

cases: u^f, = vm 3 = 1 with vm 2 = % 3 
vm 2 = v m 3 = 1-15 (all the units are arbitrary), 
the first case, the probability distribution is unimodal, 
whereas in the second case the transition to a bimodal 
distribution is observed. Indeed, the system has a phase 
transition at vm 2 — vm 3 —1.1 that corresponds to a bi- 
furcation of the critical point defined by the condition 



(15). 



In figure [2] we distinguish two regions: a drift domi- 
nated region and a diffusion dominated region. In the 
first region the chemical reactions mainly follow the gra- 
dient of the potential V(ni,n 2 ) and tend to concentrate 
around the stable critical points, so that the dynamic is 
well described by a Liouville equation 5j. In the second 
region the drift field (JTlJ is small and the fluctuations due 
to the finite size introduce a diffusive behaviour. Then 
the distribution can be approximated by the solution of 
a Fokker-Planck equation [18]. As discussed in the Ap- 
pendix, the diffusion region is approximately determined 
by the conditions 



gi{n 1 ,n 2 ) - ri(rii + l,n 2 ) 



g2{ni,n 2 ) 
0(l/N T ) 



f 1) 
(16) 



To illustrate this phenomenon, we outline in fig. [3] the 
region where condition ( 16 1 is satisfied (i.e. the gradient 
of the potential V{ni,ri2)is close to (12)). This is the 
region where the fixed points of the MM equation are 
located, and comparison with fig. [2] shows that it defines 
the support of the stationary distribution. 

In the diffusion dominated region a molecule can un- 
dergo a transition from the dephosphorylated equilib- 
rium to the double phosphorylated one. At the station- 
ary state the transition probability from one equilibrium 
to the other can be estimated by the Fokker-Planck ap- 
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FIG. 3. In grey we show the region where components of 
the vector field (111 are ~ 1 using the parameter values of 
fig. [2] (bottom). The blue lines enclose the region where the 
first component is nearby 1, whereas the red ones enclose the 
corresponding region for the second component. 



proximation; but this does not imply that the Fokker- 
Planck equation allows us to describe the transient re- 
laxation process toward the stationary state. Indeed, 
due to the singularity of the thermodynamics limit, the 
dynamics of transient states may depend critically on fi- 
nite size effects not described by using the Fokker-Planck 
approximation[T!?] . To cope with these problems in the 
PdPC model further studies are necessary. 

When the current ([8| is not zero the CME ^ relaxes 
toward a Non-Equilibrium Stationary State (NESS) and 
the field (11) is not conservative. We shall perform a 



perturbative approach to point out the effects of currents 
on the NESS nearby the transition region to the bimodal 
regime. We consider the following decomposition for the 
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vector field 



Vi by eliminating the potential A from ( 20 1 



lnai(ni,n 2 ) = -D 1 V (n 1 ,n 2 ) + D 2 H(n 1 , n 2 ) 
\na 2 {n ll n 2 ) = -D 2 V (ni,n 2 ) - DiH(ni,n 2 ) 

ni + n 2 < N T - 1 (17) 

where the rotor potential H(ni,n 2 ) takes into account 
the irreversible rotational part. The potential H(ni,n 2 ) 
can be recursively computed using the discrete Laplace 
equation 

(D 1 D 1 +D 2 D 2 )H = D 2 hx(ai(ni,n 2 ))-Diln(a 2 (ni,n 2 )) 

where n\ + n 2 < Nt — 2, with the boundary condi- 
tions H(n,N - n) = H(n,N - 1 - n) = (see Ap- 
pendix). Assuming that the potential H is small with 
respect to V, we can approximate the NESS by using 
the Maxwell-Boltzmann distribution (14) with V = Vq. 



However, as we shall show in the next section, at the 
phase transition even the effect of small currents be- 
comes critical, and the study of higher perturbative or- 
ders is necessary. To point out the relation among 
the rotor potential H, the NESS p s and the chemical 
flux J, we compute the first perturbative order letting 
Ps{n\,n 2 ) = exp(-Vo(rti,n 2 ) - Vi(ni,n 2 )). From defini- 
tion ^ the condition (8) reads: 

exp(-Vb("-i,"-2) - Vi(ni,n 2 )) n(ni,n 2 ) 

■ (exp (D 2 H(n 1 - l,n 2 )) - exp(-i3iFi(m - l,n 2 ))) 

= D 2 a(m,n 2 )) 

exp(-V r (n 1 ,n 2 ) - Vi(ni, n 2 )) r 2 {n x , n 2 ) 

• (exp (-L>iH(ni, ?i 2 - 1)) - exp (-D 2 V 1 (n 1 , n 2 - 1))) 

= -L>ia(ni,n 2 ) 

(19) 

Vi(ni,n 2 ) turns out to be an effective potential that 
simulates the current's effect on the unperturbed station- 
ary distribution by using a conservative force. We note 
that if the rotor potential H is zero, then both the cur- 
rent potential A and the potential correction V\ are zero, 
so that all these quantities are of the same perturbative 
order, and the first perturbative order of eqs. ( f 9 1 reads 



exp(-Vb(ni,n 2 ))ri(ni,n 2 ) • 

(D 1 V 1 (n 1 - l.na) + D 2 H(m - l,n 2 )) = D 2 A( ni ,n 2 ) 
exp (-Vb(ni,ra 2 )) r 2 (n\,n 2 ) ■ 

(£> 2 7i(ni,n 2 - 1) - D 1 H{n u n 2 - I)) = —DxA(ni,n 2 ) 

(20) 

From the previous equations we see that the currents 
depend both on the rotor potential H and the potential 
correction V\, which is unknown; thus they cannot be 



Di (exp(-Vo(ni,n 2 ))ri(ni,n 2 )Z>iVi(ni - l,n 2 )) 
+D 2 (exp(-Vo(^i,^ 2 ))r 2 (ni,n 2 )£) 2 Vi(ni,n 2 - 1)) 
D 2 (exp(-V (ni,n 2 ))r 2 (ni,n 2 )L»ii7(ni,n 2 - 1)) 
-L>! (exp(-V" (^i,^ 2 ))ri(ni,n 2 )D 2 lf(ni - l,n 2 )) 



(21) 



Eq. (2f ) is defined for n\ > 1, n 2 > 1 and ni + n 2 < 
Nt — 1, and we can solve the system by introducing the 
boundary conditions Vi(n,0) = Vi(0 7 n) = Vi(n,NT — 
n) = 0. It is interesting to analyze equation |2l|) in the 



phase transition regime. When the recombination terms 
ri,r 2 are almost equal and their variation is small (for our 
parameter choice this is true for ri\ ~ n 2 and ni+n 2 3> 1) 
the r.h.s. can be approximated by 

ex.p(-V (ni,n 2 ) [D 2 V Q (ni,n 2 )DxH(ni,n 2 - I) 
-DiV (n ll n 2 )D 2 H(ni - l,n 2 )] 

(22) 

As a consequence, in the transition regime this term is 
negligible since both £> 1 V r (n 1 , n 2 ) and D 2 V (ni, n 2 ) tend 
toward zero in the diffusion region where the bifurcation 
occurs; thus the first perturbative order is not enough 
to compute the stationary distribution correction, but 
higher orders should be considered. 



IV. NUMERICAL SIMULATIONS 

In order to study the non equilibrium stationary con- 
ditions in the double phosphorylation cycle we have 
perturbed the detailed balance conditions considered in 
figure (|2) by changing the value of the MM constants 
Km 2 and Km s - In figure Q we show the rotor field of 
the potential H in the cases Km 2 — Km 3 — 2.2 and 
Km 2 = Km 3 = 1-8 when the detailed balance condition 
( 13 1 does not hold. We see that in the first case (left 



picture) the rotor field tends to move the particles from 
the borders toward the central region n\ — n 2 , so that 
we expect an increase of the probability distribution in 
the center, whereas in the second case the rotor field is 
directed from the central region to the borders and we 
expect a decrease of the probability distribution in this 
region. To illustrate the effect of the potential H, we 
compare the zero-order approximation of the probabil- 
ity distribution ( |f 4[ ), where the potential Vo(ni,n 2 ) is 
computed using decomposition (17) with the stationary 
solution of the CME The main parameter values are 
reported in table [TV| 

Case Km 2 K Ms vm 2 vm 3 

I 1.8 1.8 1.05 1.05 

II 1.8 1.8 1.15 1.15 

III 2.2 2.2 1.05 1.05 

IV 2.2 2.2 1.15 1.15 



directly computed from ( 20 ) . We obtain an equation for 
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FIG. 4. Plot of the rotor field for the potential H using the following parameter values: vmi = 1>M2 = 1, vm 2 = u a/ 3 = 1.15 
Km-i = Km 4 = 1, Nt = 40 and Km 2 = Km 3 = 1.8 (left picture) or Km 2 = Km 3 = 2.2 (right picture). 




FIG. 5. (Left picture): plot of the zero-order approximation for the probability distribution using the decomposition (17 1 for 
the vector field associated to the CME. (Right picture): plot of the stationary distribution computed by directly solving the 
CME p|. We use the parameter values of the case I in the table IIVI 



whereas the other parameters values are: Vm x = % 2 = 
1,Km! — Km 4 — 1 and Nt = 40. For the first parameter 
set, the zero-order approximation is a bimodal distribu- 
tion, but the effect of the currents induced by the rotor 
potential H (see fig. [4]) left) are able to destroy the bi- 
modal behaviour by depressing the two maximal at the 
border (see fig. 5). 

If we increase the value of the K^m and K$m (case II) , 
the exact stationary distribution becomes bimodal and 
the effect of currents is to introduce a strong transition 
probability between the two distribution maxima (fig. [6|. 

Therefore, when the MM constants Km2 and Km3 are 
< 2 (we note that for Km 2 = Km 3 = 2, the detailed 
balance holds), the non-conservative nature of the field 
(17) introduces a delay in the phase transition from a 



mono-modal to a bimodal distribution. However, when 
we consider Kmi = Km3 > 2 (cases III and IV) the rotor 
potential H moves the particle towards the borders and 
the central part of the distribution is depressed. This is 
shown in the figure [7] where we compare the zero-order 
approximation of the stationary distribution and the so- 
lution of the CME using the case III parameters of the 
table HVl 



Finally, in case IV of table |IV[ the CME stationary 
solution undergoes a transition to a bimodal distribu- 
tion, whereas the zero-order approximation is still mono- 
modal ( fig. [9|, so that the effect of currents is to antici- 
pate the phase transition. 

Using the stationary solution one can also compute the 
currents according to definition ([5| . In figure (9) we plot 



FIG. 6. The same as in fig. [5] using parameter values of case II in table |IV| 
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FIG. 7. The same as in fig. [5] using the parameter values of the case III in table |TV] 



the current vector in case IV parameters to show that 
the current tends to become normal to the distribution 
gradient near the maximal value. 

This result can be also understood using the perturba- 
tive approach ( 21 ), where one shows that the main effect 



of the V\ potential correction is to compensate the rotor 
field of H along the distribution gradient directions. As a 
consequence, the current is zero at the maximal distribu- 
tion value and condition (151 defines the critical points of 



the stationary distribution even in the non-conservative 
case. 



V. CONCLUSIONS 

The CME approach we present here is a powerful 
method for studying complex cellular processes, even 
with significant simplifications such as spatial homogene- 
ity of volumes where the chemical reactions are taking 
place. The CME theory is attractive for a variety of 
reasons, including the richness of aspects (the capabil- 



ity of coping with fluctuations and chemical fluxes) and 
the possibility of developing thermodynamics, starting 
from the distribution function. The violation of detailed 
balance gives information on the "openness" of the sys- 
tem and on the nature of the bistable regimes, which are 
induced by the external environment; in contrast, it is 
a free-energy equilibrium when detailed balance holds. 
This statement can be expressed in a more rigorous form 
by introducing the vector field generated from the ra- 
tio between the generation and recombination terms, by 
decomposing it into a sum of "conservative" and "rota- 
tional" fields (Helmholtz decomposition) and by relating 
the chemical fluxes to the non-conservative field. The 
magnitude of deviations from detailed balance influences 
the form of the stationary distribution at the transition to 
a bistable regime, which may be driven by the currents. 
An interesting test for the prediction of the PdPC CME 
model would be to perform one experiment with the pa- 
rameter values chosen to satisfy DB and compare it with 
another set of parameters where DB is not fulfilled. Our 
results show that the PdPC can operate across these two 
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FIG. 8. Current vector field computed by using definiton 
|5| and the stationary solution of the CME with case IV pa- 
rameters. The distribution is bimodal (cfr. figure 8) and the 
current lines tend to be orthogonal to the distribution gradi- 
ent near the maximal value. 



regions, and that the transition regime can be explained 
by the role of the currents, that, within a thermody- 
namic framework, can be interpreted as the effect of an 
external energy source. A full thermodynamic analysis 
of this cycle is beyond the scope of this paper, but we can 
surmise that this approach might be extended to other 
cycles in order to quantify if, and how much energy, is 
required to maintain or create bistability. Another way 
to extend this analysis would be the generalization to 
n-step phospho/dephosphorylation cycles, where the sta- 
tionary distribution will be the product of n independent 
one-dimensional distributions. In conclusion, our results 
could be important for a deeper characterization of bio- 
chemical signaling cycles that are the molecular basis for 
complex cellular behaviors implemented as a "switch" 
between states. 
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FIG. 9. The same as in fig. [5] using parameter values of case IV in table |IV| 



VI. APPENDIX 

The master equation describes the evolution of one- 
step Markov Processes according to 



(ni,n 2 ,t) 



dp 

at 

01 (ni - l,n 2 )p(n 1 - l,n 2 ,t) - g 1 (n 1 ,n 2 )p(n ll n 2l t) 
+rt(m + l 1 n 2 )p(ni + l,n 2 ,t) - n(ni, n 2 )/?(ni, n 2) t) 
+g 2 (ni,n 2 - l)p(n u n 2 - l,i) - g 2 (ni, n 2 )p(ni, n 2) i) 
+r 2 (m,n 2 + l)p(ni,n 2 + l,t) - r 2 (ni, n 2 )/?(ni, n 2) i) 

(23) 

with the boundary conditions for the coefficients 



gi(n, N — n) = g 2 (n, N — n) = and 
r 2 (n,0) =ri(0,n) = ne[0,7V T ] 



(24) 



so that ?ii + n 2 < JVy. By introducing the difference 
operators |7]), eq. (23) can be written in the form of a 
continuity equation 



dp 
dt 



(ni,n 2 ,<) = -D 1 Ji(m,n 2 ,t) - D 2 J 2 (n ll n 2l t) (25) 



where we introduce the current vector J of components: 

Ji(m,n 2 ,t) = 5i(ni - l,n 2 )p(m - l,n 2 ,t) 
-ri(ni,n 2 )p(ni,n 2 ,i) 

J 2 {n u n 2 ,t) = g 2 {n u n 2 - l)p(m,n 2 - l,i) 
-r 2 (ni,n 2 )p(ni,n 2 ,i) 

(26) 

The stationary solution p s (ni, n 2 ) is characterized by the 
zero divergence condition for the current (25). Detailed 



balance holds when the current is zero and p s satisfies 

r 1 {n 1 ,n 2 )p s (n 1 - l,n 2 ) 

Ps(ni,n 2 ) gx(nx - l,n 2 )\ _ () 



= 



p s (ni- l,n 2 ) n(ni,n 2 ) 
r 2 {ni,n 2 )p s (ni,n 2 - 1) 

Ps{ni,n 2 ) _ g 2 (ni,n 2 - 1) 
p s (n 1) n 2 - 1) r 2 {ni,n 2 ) 



(27) 

for < ni and < n 2 . The previous equations can be 
written in the form 



L>iln(p s (ni - l,n 2 )) = In 
D 2 \n(p s (m,n 2 - 1)) = In 



gi(ni - l,n 2 ) 
ri(ni,n 2 ) 

g 2 {ni,n 2 - 1) 
r 2 (n 1 ,n 2 ) 



(28) 



and, if one introduces the the vector field 

£i(ni,n 2 ) 



ai(ni,n 2 ) 
a 2 (ni,7i 2 ) 



n(ni + l,n 2 ) 

92(ni,n 2 ) 
r 2 (ni,n 2 + 1) 



(29) 

due to the commutative property of the difference opera- 
tors Di , detailed balance implies an irrotational character 
for the vector field ln(a(ni, n 2 )) 



D 2 ln(ai(m, n 2 )) - D 1 ln(a 2 (ni, n 2 )) = 



(30) 



If we have no singularities in the domain, eq. ( 30 ) is a suf- 



ficient condition for the existence of a potential V{ni, n 2 ) 
(cfr. eq, (12)) and the distribution p s (ni, n 2 ) can be com- 



puted using the recurrence relations 

p s {ni + l,n 2 ) = ai(ni,n 2 )p s (ni,n 2 ) 
p s (ni,n 2 + 1) = a 2 (ni,n 2 )p s (ni,n 2 ) 



(31) 
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Therefore, the components ai (711,712) and 02(711,712) can 
also be interpreted as creation operators according to 
relations plj) and detailed balance condition pOl is 



p s (n\, 71%) by applying recursively the relations (31) in 



equivalent to the commutativity property for these op- 
erators. The stationary distribution can be written 
in the Maxwell-Boltzmann form ( 14 ) and the potential 
V (711,712) is associated with an "energy function". We 
finally observe that the critical points of the stationary 
distribution p s are defined by the condition 



01(711,712) = a 2 (ni,n 2 ) = 1 



(32) 



For the double phosphorylation cycle (JlJ it is possible to 
derive explicit expressions for the stationary distribution 



a specific order: for example, first moving along the n 2 
direction and then along n\, we obtain an expression for 
p s (ni,7i2) as a function of p s (0,0) 



I \ TT 9i{i ~ 1>«2) tt 92(0,1 - 1) 

ps{ni > n2) = n ri(i ,„ 2) n w ft(0, 0) 

(33) 

A direct substitution of the coefficients Q in the relation 
(33) gives 



/ \ tt K M4 V M 2(N T - i - n 2 + 1) K M1 K M3 + K M1 (N T - i - n 2 ) + K M3 i 

Ps (7li,7l2) = I I -. r ' 

xl . K M 2K M 4 + K M a(N t - 1 - n 2 + 1) + K M 2n 2 K M3 Vm^ 

tt K M1 V M3 (N T -l + l) K M2 KM4 + K M i(N T -l)+KM2l { ) 
'\1KmiK M3 + Kmi(N t -1 + 1) K M2 V MA l Ps[ ' 



r 



We can further simplify this expression by using the 
definition of multinomial coefficients and the rising and 
falling factorial symbols, defined as x^ 71 ' = x(x + l)(x + 



2)- 
2)- 



(x 
(x 



■1) 
1) 



(x+n—l)\ 
(x-1)! 



and X( n ) = x(x — l)(x 



(x—n)\ 



respectively. 



Ps(ni,n 2 ) 

Km3 



( Vm3 



( Nt — n 2 \ ( Nt 



Vmi J \Vm4 
Kmi \ ni ( Km2 — K 



M4 



\ 7ii J \ n 2 

Tl2 



Km3 J \ K M 2 

(K M1 (l + n 2 - K M3 ~ N T ) - KM3) {ni) 
(K M i ~ K M3 )M(K M2 {1 + ^) + N T - n 2 ) (ni) 
(K M2 (\ + Kma) + K Mi (N T - !))("=) 



(K M 2 - K M4 )(^)(Km3 + N T )(n 2 ) 



Finally, it is interesting to go to a continuous limit 
that is equivalent to Nt — > 00. First we introduce the 
population densities A — ni/Nx and A p = ti 2 /Nt and 
use the fact that the generation and recombination rates 
are invariant by substituting n\ and n 2 with A and Ap. 



Then we approximate 

D 1 V(A, Ap) = V(A- 



1 dV 
NrdA 



(A 



l/N T ,A P p)-V(A,A p P ) 
f~^) + 0(l/iV|) 



and a similar expression holds for D 2 V(A, Ap). Accord- 
ing to eq. (12), the partial derivatives of V(A, Ap) are 
bounded when Nt —> 00 only in the domain where the 
following approximation holds (diffusion dominated re- 
gion) 



g l (AAp") 
n(A,A p ) 



1 + 0(1/N 7 



1,2 



(34) 



and we can estimate 



In 



9i(AA p ) 



n(A + l/N T ,A p ) 
t n(A + l/N T ,A$)- 



9i(A,A p ) 



In 



r x {A+l/N T ,A p )+g x (A,A p ) 
92(A,A p ) 



0(l/N 



r 2 (A,A p + l/N T ) 
r 2 (A,A p + l/N T )- 



92(A,A p ) 



J r 2 (A,A p + l/N T )+g 2 (A,A p ) 



0(l/N 



(35) 



Then we may approximate (we use the convention of leav- 
ing out the dependence on Ap) 
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r x (A + 1/N T ) - 51 (A) _ n(A + 1/2N T ) - gi(A + l/2N T ) + 1/2N T {dri/dA[A + l/2N T ) + d gi /dA(A + l/2N T )) 
n(A+l/N T )+ gi (A) ~ ri (A + l/N T )+ gi {A) 

(36) 



up to an error of order 0(1/N T ) (a similar expression is 
obtained for the second equation). Therefore, detailed 
balance in the continuous limit reads 



n(A, A P P ) - gi (A, A P P ) + l/2N T (dr 1 /dA(A,A p )+dg 1 /dA(A,A p )) ^ 1 dV P 
r 1 (A + l/2NT,A p )+g 1 (A-l/2N T ,A p ) 2N T 3A [ ' p) 

r 2 (A,A p )-g 2 (A,A p ) + l/2N T (dr 2 /dA p (A,A p )+dg 2 /dA p (A,A p )) _ 1 dV 



r 2 (A, A p + 1/2N T ) + g 2 (A, A p - l/2JV r ) 2N T dA p 



{A, Ap) 



(37) 



I 

The limit Nt — > oo turns out to be singular since 



-?Y-(A,A P ) 



2N T ( ri (A, A p ) - gi (A, A p )) + dn/dA(A, A p ) + d 9l /dA(A, A p ) 
ri (A,A p ) + gi (A,A p ) 



dV 

oaJ 



, a p 2N T (r 2 (A, A p ) - g 2 (A, A p )) + dr 2 /8A p (A, A p ) + dg 2 /dA p (A, A p ) 
[ ' P> r 2 (A,A p )+g 2 (A,A p ) 



(38) 



Hence in the diffusion domain defined by condition ( 34 1 , fields 



we recover detailed balance for a Fokker-Planck equation 
with drift and diffusion coefficients defined as: 

d(A, Ap) = 2N T ( n (A, Ap) - 9i (A, Ap)) i = 1,2 



and 



b i (A,A p ) = r i (A,A p )+g i (A,A p ) 



i = 1,2 



In the diffusion region the drift and the diffusion coef- 
ficients are of the same order, otherwise Ci(A,A p ) 3> 
bi(A,Ap) when Nt ^> 1. As a consequence, the sta- 
tionary solution of F.P. equation is an approximation of 
the stationary distribution of the CME in the diffusion 
region, but the approximation of the transient state dy- 
namics using the F.P. equation requires further studies 
due to the singularity of the thermodynamic limit. 
In the generic case of eq. (23 1, we represent the r.h.s. of 
cq. ( 28 1 as a sum of a rotational and a gradient vector 



lnai(ni,n 2 ) = DiV(ni,n 2 ) + D 2 H (m, n 2 ) 
\na 2 (ni,n 2 ) = D 2 V{ni,n 2 ) - DiH(m,n 2 ) 

Taking into account the condition ni+n 2 < Nt — 1, from 
the eqs. (39) we get the discrete Poisson equations 



(D 1 D 1 +D 2 D 2 )V(n 1 ,n 2 ) = 
Di ln(ai(m,n 2 )) + D 2 In (a 2 (m,n 2 )) 

{D 1 D 1 +D 2 D 2 )H(n 1 ,n 2 ) = 
D 2 In (ai(ni,n 2 )) - D\ In (a 2 (m,n 2 )) 



(40) 



We remark that the r.h.s. of eqs. (40) is defined only if 
rii + n 2 < N — 2 and corresponds to N(N — l)/2 inde- 
pendent equations, whereas we have (N + 2) (TV + l)/2 
unknown values H(ni,n 2 ). As a consequence from the 
explicit form of the discrete Poisson operator 



(Dj + D 2 2 )H{ nil n 2 ) = H(m + 2, n 2 ) ~ 2H(n x + 1, n 2 ) + H(m,n 2 ) + H{n u n 2 + 2) - 2H{n u n 2 + 1) + H{m,n 2 ) 

I 

we can set the boundary conditions H (n, N — n) — setting 

H(n, N — 1 — n) = and recursively solve the system 2H(n N — 2 — n) = 

Oi(n, N - 1 - n)a 2 (n, N-2-n) 



In 



a 1 (n,N-2 - n)a 2 (n + l,N- 2-n) 



(41) 
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and successively using the equations 



2H(n, N-j-n) = -H(n + 2,n-j-n)- H(n, N-j-n-2) + 2H(n + 1, N - j - n) + 2H(n, N-j-n-l) 

. ai(n, N - j -n + l)a 2 (n,N - j - n) 
ai(n, N — j — n)a 2 (n + 1, N — j — n) 

(42) 



for N > j > 2. Once H{n\,n 2 ) is computed, we define 
the "potential" V {111,112) by using eq. (39). The recur- 



sion relations (|42j) can be written in an exponential form 
by defining 



R(m, n 2 ) = exp(i?(ni, n 2 )) 



R(n,N-j-n)) 



Ux{n,N- 


- j 


-n+l)a 2 (n,iV 


-j 


-n) 


V oi(n,JV- 




-n)a 2 (n+l,iV 


-j 


-n) 


ii(n + l,JV- 


-j - 


- n)iZ(n, N — j - 


- 11 - 


-1) 


+ 2,?i 


- J 


- n)R(n,N - j 


— n 


-2) 



(43) 



As a consequence, the recurrence (31) reads 



p s {ni + l,n 2 ) = ai(ni,n 2 )— — — p s (n 1 ,n 2 ) 

R{n 1 ,n 2 + 1) 

p s {ni,n 2 + 1) = a 2 (ni,n 2 ) — — r— p s (n 1 ,n 2 ) 



R(ni,n 2 ) 



(44) 



for all rii +n 2 < JV - 1. 

The current componentss (25) turn out to be propor- 
tional to the rotational part of the field (39) (i.e. to 



H(ni . no)) [11]. so that the current vanishes at the points 
where condition (32) is satisfied. One can prove that 



the critical points of the stationary distribution are still 
defined by eq. (32). Indeed, if one computes the for- 



mal expansion of the generation and recombination rates 



around a solution of eqs. (31) 

9g, 



g(n) = l + ^(n*)-An + .. 
dr 

r{n) = 1 + g^(n*) ■ An + .. 



(45) 



(for the sake of simplicity we have normalized the value of 
the generation and recombination rate to 1 at the critical 
point) the current components can be approximated by 
the expressions 

Ji(ni,n 2 ) ~ p s (m - l,n 2 ) - p s (nx,n 2 ) 



+ -^-( n *) ■ Anp s (m - 1,712) .. 
on on 



(n*) ■ Anp s (ni,n 2 ) 



J£(ni,n 2 ) ~ p s {ni,n 2 - 1) - p s (n 1 ,n 2 ) 



9.92 



dr 2 
On 



(n*) • Anp s {n l7 n 2 ) 
(46) 



At the critical point n* we get 

^(ra* - 1,713) = Ps{nl,n* 2 - 1) = p s (ni,7i 2 ) (47) 



since J s (71* , 712 ) = 0. The condition ( |47| ) means the n* is 
a critical point for the stationary distribution p s . 
When H(ni,n 2 ) is small, the detailed balance solution 
( 14 ) is a good approximation of the stationary solution 



p~Jrii,n 2 ) and a perturbative approach can be applied. 
Let us write the stationary condition Q in the form 

Psjni - l,n 2 ) 

ps(ni,n 2 ) 
p s (n 1 ,n 2 - 1) 
p a (711,712) 



Diri(m,n2)Ps(ni,n 2 ) yL - oi(ni - l,n 2 ) 
D 2 r 2 (m,n 2 )p s (ni,n 2 ) 1 - a 2 (ni,n 2 - 1) 



By using the definitions (39), we assume that the rota- 



tional field is associated with a potential eH(ni, n 2 ), with 
e«l perturbation parameter and we write the station- 
ary solution in the form 

Ps{ni,n 2 ) = Cexp(V(n 1) n 2 ) + eV 1 (n 1 , n 2 )) (49) 

From a direct calculation we get 
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D 1 r 1 (n 1 ,n 2 )e v ^^'> (1 - exp(e(£> 2 # (m - l,n 2 ) - £>i^(ni - l,n 2 )))) 
+D 2 r 2 (n 1 , n 2 )e v{n ^ (1 - exp(-(e£> 1 fl'(ni, n 2 - 1) + £> 2 Fi(ni,n 2 - 1))) ~ 
e J D 1 r 1 (n 1 ,n 2 )e v '(™ 1 ^ 2 ) (DiV^m - l,n 2 ) - D 2 H(m - l,n 2 )) 
+e^ 2 r 2 (n 1 ,n 2 )e y( ™ 1 '" 2) (£> 2 Vi(»i, n 2 - 1) + £>iif(m, n 2 - 1)) = 

(50) 

I 

for all the values n\ + n 2 < N — 1 and > 1. The the previous equation, and it enters in the definition of 
correction potential V\{ni,n 2 ) has to be computed from the stationary currents. 



